Dynamical Power Spectra (on real data)


In [1]:
%matplotlib inline

In [2]:
# load auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits

# import stingray
import stingray

plt.style.use('seaborn-talk')

All starts with a lightcurve..

Open the event file with astropy.io.fits


In [3]:
f = fits.open('emr_cleaned.fits')

The time resolution is stored in the header of the first extension under the Keyword TIMEDEL


In [4]:
dt = f[1].header['TIMEDEL']

The collumn TIME of the first extension stores the time of each event


In [5]:
toa = f[1].data['Time']

Let's create a Lightcurve from the Events time of arrival witha a given time resolution


In [6]:
lc = stingray.Lightcurve.make_lightcurve(toa=toa, dt=dt)

In [7]:
lc.plot()



DynamicPowerspectrum

Let's create a dynamic powerspectrum with the a segment size of 16s and the powers with a "leahy" normalization


In [8]:
dynspec = stingray.DynamicalPowerspectrum(lc=lc, segment_size=16, norm='leahy')

The dyn_ps attribute stores the power matrix, each column corresponds to the powerspectrum of each segment of the light curve


In [9]:
dynspec.dyn_ps


Out[9]:
array([[  2.01901704e+00,   2.32485459e+00,   5.14704363e+00, ...,
          9.76872866e-01,   9.49269045e-01,   4.60522187e+02],
       [  2.93960257e+00,   2.48892516e+00,   3.39280288e+00, ...,
          6.23511732e+00,   4.27550837e+00,   1.06261843e+02],
       [  3.64619904e+00,   1.58266627e+00,   3.42614944e-01, ...,
          1.16952148e+00,   3.54994270e+00,   4.56956463e+01],
       ..., 
       [  1.69311108e+00,   5.18784072e-01,   1.57151667e+00, ...,
          1.09923562e+00,   3.40274378e-01,   2.53108287e+00],
       [  2.95675687e-01,   2.47939959e+00,   2.84930818e+00, ...,
          2.99674579e-01,   1.48585951e+00,   7.49068264e+00],
       [  8.84156884e-01,   1.65514790e+00,   4.17385519e-01, ...,
          7.54942692e+00,   9.99801389e-01,   2.03835451e-01]])

To plot the DynamicalPowerspectrum matrix, we use the attributes time and freq to set the extend of the image axis. have a look at the documentation of matplotlib's imshow().


In [10]:
extent = min(dynspec.time), max(dynspec.time), max(dynspec.freq), min(dynspec.freq)

plt.imshow(dynspec.dyn_ps, origin="lower left", aspect="auto", vmin=1.98, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)


Out[10]:
(700, 850)

In [11]:
print("The dynamical powerspectrun has {} frequency bins and {} time bins".format(len(dynspec.freq), len(dynspec.time)))


The dynamical powerspectrun has 65535 frequency bins and 104 time bins

Rebinning in Frequency


In [12]:
print("The current frequency resolution is {}".format(dynspec.df))


The current frequency resolution is 0.0625

Let's rebin to a frequency resolution of 2 Hz and using the average of the power


In [13]:
dynspec.rebin_frequency(df_new=2.0, method="average")

In [14]:
print("The new frequency resolution is {}".format(dynspec.df))


The new frequency resolution is 2.0

Let's see how the Dynamical Powerspectrum looks now


In [15]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=1.98, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(500, 1000)


Out[15]:
(500, 1000)

In [16]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)


Out[16]:
(700, 850)

Rebin time

Let's try to improve the visualization by rebinnin our matrix in the time axis


In [17]:
print("The current time resolution is {}".format(dynspec.dt))


The current time resolution is 16.0

Let's rebin to a time resolution of 64 s


In [18]:
dynspec.rebin_time(dt_new=64.0, method="average")

In [19]:
print("The new time resolution is {}".format(dynspec.dt))


The new time resolution is 64.0

In [20]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)


Out[20]:
(700, 850)

Trace maximun

Let's use the method trace_maximum() to find the index of the maximum on each powerspectrum in a certain frequency range. For example, between 755 and 782Hz)


In [21]:
tracing = dynspec.trace_maximum(min_freq=755, max_freq=782)

This is how the trace function looks like


In [22]:
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', alpha=1)
plt.show()


Let's plot it on top of the dynamic spectrum


In [23]:
extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
           interpolation="none", extent=extent, alpha=0.7)
plt.colorbar()
plt.ylim(740,800)
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', lw=3, alpha=1)
plt.show()


The spike at 400 Hz is probably a statistical fluctutations, tracing by the maximum power can be dangerous!

We will implement better methods in the future, stay tunned ;)